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We consider growing spheres seeded by random injection in time and space. Growth stops when 
two spheres meet leading eventually to a jammed state. We study the statistics of growth limited by 
packing theoretically in d dimensions and via simulation in d = 2, 3, and 4. We show how a broad 
class of such models exhibit distributions of sphere radii with a universal exponent. We construct a 
scaling theory that relates the fractal structure of these models to the decay of their pore space, a 
theory that we confirm via numerical simulations. The scaling theory also predicts an upper bound 
for the universal exponent and is in exact agreement with numerical results for d — 4. 
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I. INTRODUCTION 

The dynamic packing of objects is an often overlooked 
variation on the theme of static packing. Given a mech- 
anism for the creation, growth, movement, and interac- 
tion of like objects in a given dimension, what struc- 
tures result? Here, wc find universal features of a sim- 
ple yet broad class of such mechanisms falling under the 
rubric "packing-limited growth" (PLG). A PLG mech- 
anism entails objects being seeded randomly, growing 
according to a rule that may be specific to each object, 
and stopping when a part of their boundary hits that of 
another object. A motivating physical example of this 
kind of pattern formation may be found in the compe- 
tition between tree crowns in dense forests [0, |||, the 
structure of porous media |j, |], and the generalized 
problem of dense packings [p|. 

A simple model of PLG has previously been studied by 
Andrienko, Brilliantov, and Krapivsky 0, ^ (the ABK 
model). In their setting, spheres are seeded randomly in 
space and time. A sphere's radius increases linearly with 
time, halting when another sphere is touched. This mod- 
el is amenable to an approximate analysis for d > 1 || 
and has an exact solution in d = 1. In this present work, 
we are interested in the limiting distribution of radii, 
N(r). For PLG models, we expect N(r) oc r~ a for small 
r. Note that the fractal dimension D of the set compris- 
ing all sphere centers, which is often measured instead of 
a, is related to a as 13 = a — 1 B. In d = 1, the exact 
solution gives a = 1, which corresponds to D = (mean- 
ing the number of centers diverges logarithmically) ||. 
For d — 2, Andrienko et al. j?J report that D ~ 1.75 
and hence a ~ 2.75 based on numerical evidence. Else- 
where, the same authors §| determine numerically that 
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a ~ 2.53. 

We claim that the actual value of a is essentially inde- 
pendent of the specifics of the growth dynamics. To see 
why, we examine a related d — 2 random packing mod- 
el due to Manna ||. The packing process begins in a 
finite-sized volume (or alternatively with an initial pop- 
ulation of fixed radius disks). New disks are added one 
at a time by randomly choosing a point in the packing's 
pore space and centering there the largest possible non- 
overlapping disk. Manna finds a ~ 2.62 and a ~ 2.64 for 
two example packings. 

We refer to the Manna model as "random Apollonian 
packing" (RAP) since it may be seen as a variation of 
the well known Apollonian packing For the d = 2 
version of the latter, pore spaces are always formed by 
three disks each touching the other two and disks are 
added so as to fill the pore space fully (i.e., the insert- 
ed disk touches all three of its surrounding neighbors). 
Numerics give a ~ 2.31 for Apollonian packing in agree- 
ment with analytically derived bounds [[IT] . Aste (T^] has 
shown for static polydisperse packings that are space fill- 
ing, Apollonian packing provides the lowest bound on a 
while the upper bound is a = d + 1 . 

There are two important observations to make here. 
First, as noted by Brilliantov et al. J8|, the RAP model is 
the ABK model in the case of infinitely fast growth: as 
soon as a sphere is nucleated, it instantaneously expands 
to hit the nearest sphere boundary. Second, for general 
PLG models, pore spaces evidently increase in number 
and decrease in size with time. This means that the like- 
lihood of two spheres nucleating in the same pore space 
also decreases. Eventually the mechanism of the RAP 
model must take over, and all collisions will be between 
a growing sphere in a pore space and a "stuck" sphere 
forming part of the pore space boundary. Thus, the RAP 
model is not just a curious end point of the ABK model 
but, in fact, entails the sole mechanism describing how 
small radius spheres pack. We suggest then that mea- 
surements of a in the ABK model should coincide with 
revised estimates from simulations of the RAP model. In 
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the remainder of the paper, we describe a general class 
of PLG models for which a is invariant, provide a simple 
solution for the d = 1 problem, develop a scaling theo- 
ry that describes how the radii distribution and volume 
of pore space evolve with time, and provide results from 
extensive numerical simulations. 



II. GENERAL MODEL OF PACKING-LIMITED 
GROWTH 



For our general conception of PLG we take d- 
dimensional spheres growing in a volume V. Spheres are 
nucleated randomly in space at a rate k per unit volume 
surviving only if they are injected into the unoccupied 
pore space. All spheres stop growth upon contact with 
a neighboring sphere. Figure [j] provides a visual context 
for the dynamics explained below. The ith sphere, which 
is initiated at ti, grows at a rate Gi{t — U) giving the 
radius ri(t) as 

rt— u 

n(t) = / Gi(u) du, (1) 
Jo 

for t > ti. We assume the weak requirement that each 
sphere grows in a strict monotonic fashion, i.e., for each 
i, Gi{t — ti) > € > 0. The model is thus one of arbitrary 
individual growth limited solely by packing. The dynam- 
ics continue in the above fashion until the volume V is 
completely filled and a final jammed state is reached. 

Consider an individual pore T of diameter A (i.e., the 
maximum separation of two boundary points that can be 
joined by a line contained entirely within the pore). The 
rate v at which spheres nucleate in T (excluding growth 
for the moment) is given by 

v s KV d X d , (2) 

where Vd is the volume of <i-dimensional sphere of unit 
radius. The typical time r between nucleations in T is 
therefore 

r = 1/v ~ K-*V d - x X- d . (3) 

On the other hand, if sphere i nucleates in F, it will 
jam in a time rj am that can be no greater than the time 
it takes for its radius to reach A/2. Together with the 
assumption that Gi(t — ti) > e > 0, this gives 

Tjam < A/2e. (4) 

So when rj am <C r, i.e., when 

V d X d+1 K/2e < 1, (5) 

we expect the packing mechanism to be the same as the 
RAP model — all spheres will be stopped by an exist- 
ing stopped sphere and never by another moving sphere. 
Growth rate therefore becomes irrelevant as far as the 
final packing is concerned. From Eq. (@), we see that 



(c) (d) 

FIG. 1: The dynamics and limiting states of PLG models. 
In (a) and (b), collisions between centers are marked with 
lines while in (c) and (d), the region contained by spheres is 
black and the pore space is left white. All four pictures come 
from the same simulation, (a) and (c) depict the condition 
of a subsection of a system packed with 500 spheres, (b) and 
(d) depict the same subsection after 2000 spheres have been 
packed. 



the general model reduces to the RAP model when the 
maximum pore size A max satisfies 

A max « (2e/V d K) 1/(d+1) . (6) 

Since A max steadily decreases with t (i.e., space fills up) 
this condition will always be eventually satisfied. 



III. EXACT d = 1 SOLUTION 

As we have noted, the ABK model has been solved 
exactly in the d = 1 case || with the outcome being 
a = 1 . Here, we achieve the same result with a consider- 
ably simpler calculation. Since we have posited that a is 
a universal exponent for a general class of packing-limited 
growth models, we may choose the straightforward exam- 
ple of the d = 1 RAP model. We take the unit interval 
[0, 1] as the initial vacant pore space. Spheres are now 
line segments and are limited either on the left or right as 
they are added with the points and 1 providing initial 
boundaries. Thus, the unit interval is filled in from each 
end with the one solitary pore space diminishing in the 
middle. Writing the length of the nth line segment as l n 
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we have 



n-1 



In Z n —i I 1 ^ 1% J ; 



where Zi is a random number uniformly distributed on 
the unit interval and l\ = zq. In other words, each new 
line segment is a random fraction of the current pore 
space. Iteratively solving Eq. (Q) gives 



(8) 



i=0 



We can in principle find the distributions of the li but 
for our present purposes their means are sufficient and 
we have 



(9) 



We therefore expect the typical number of line segments 
with I' > I = 2~ n to be 



tf(l'>l = 2-)=» = -£i. 



(10) 



Since this is the cumulative frequency distribution, i.e., 
N(l' >l) = J™ N(l')dl', we find the frequency distribu- 
tion N(l) must behave as 




FIG. 2: P(r;n) for a d = 2 RAP simulation with n = 10 6 . 
The inset distribution is Pins(r; which is obtained by simu- 
lating the RAP procedure for another 10 6 disks without actu- 
ally inserting any of them into the pore space. The dashed 
line in the inset figure corresponds to the theoretical predic- 
tion P ins (0;n) = S(n)/$(n). 



N(l) * -r\ 



(11) 



yielding, as expected, a = 1. 



IV. DESCRIPTION OF INSERTION 
PROBABILITY AND PORE SPACE VOLUME 

We next investigate the form of Pi ns (r;n), the prob- 
ability distribution of the (n + l)th sphere's radius to 
be inserted into a packing. In addition, we characterize 
P(r; n), the probability distribution of sphere radii after 
n spheres have been packed. Note that since we have 
shown that a general class of PLG models reduce to the 
RAP model, we now use the number of spheres n rather 
than time t to index our quantities. 

We first note that the distribution P(r; n) fills in from 
the right with increasing n. As the packing fills in 
space, the maximum pore size either decreases or does 
not change and so does the maximum size of any sphere 
that may be added. We write the radius of this largest 
sphere as r c . We observe that to a first approximation, 
P(r; n) above this cutoff scale r c follows its limiting power 
law form while below the distribution it is essentially flat 
(see Fig. ||) . For the purposes of estimation, we assume 
this form exactly as 



P(r;n)= a ° 



for r < r c 
for r > r r 



(12) 



The corresponding frequency distribution is N(r; n) oc 
P(r;ri). Since the tail of N(r;n) remains fixed as n 
increases, we can obtain the scaling of r c with n. From 
Eq. (p~3), the tail of N(r; n) behaves as 



N(r;n) 



a — 1 



a 



= kr 



(13) 



for r > r c . Since the prefactor k must be constant, we 
have 



a-1 

— TT' 1 
ak 



-l/(a-l) 



-l/(a-l) 



(14) 



The uniformity of P(r; n) for r < r c closely relates to 
the form of Pi ns (r;n). It is possible to write down an 
exact expression for Pins(r; n), that is, 



Pins(r;n) 



J dVS(D n (x) - r) 
JdV 



(15) 



where the integrals are over the pore space, and D n (x) 
is the distance from the point x to the closest pore space 
boundary after n spheres have been inserted. This inte- 
gral may be solved exactly in the limit of very small radii, 



lim P ias (r;n) 

r—>0 



S(r 



$(n)' 



(16) 



where S(n) and $(n) are the surface area and available 
pore space of the existing n packed spheres. This means 
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that the region of pore space available for insertion of an 
infinitesimal sphere is proportional to the surface area of 
the extant packing. Assuming that this holds approxi- 
mately for all spheres below the cutoff, the full distribu- 
tion Pins(^; n) may be modeled as a purely flat distribu- 
tion, 



Pms(r;n) 



c 1 for r < r c 
for r > r r . 



(17) 



Good agreement with this approximation is shown in the 
inset of Fig. ||. Comparing Eq. jl6| ) with Eq. ( fl7| ) , we see 
that r c — $(n)/S(n). We calculate S(n) and 3>(n) using 
N{r;n), 



S(n) 



and 



kaVd 



kaV d 



-( a -d) 



y ,d+l-a 



(d + l)(d + 1 - a) 



(18) 



(19) 



where as before Vd is the volume of a unit radius sphere 
in d dimensions. Note that S(n) — > oo and $(n) — > as 

71 — > CO. 



Using Eqs. (18) and (|19[), and the result r c 
$(n)/S(n), we find an estimate of a as 



d + 2 



(20) 



In contrast, a modified calculation for the ABK model 
finds l§ 



(XmS = 1 + d ( 1 + exp 



2 d+2 _ 2 

d+2 



(21) 



Note that a must be an upper bound on the true val- 
ue of a for normal Apollonian packing. Furthermore, 
due to the nature of the approximations made in form- 
ing Pi ns (r;n) and P(r;n), a is also an upper bound for 
the exponent of the RAP model. As we show in the fol- 
lowing section when we consider our numerical results, 
both of these predictions of a appear to hold for certain 
(mutually exclusive) ranges of d. 

Using Eqs. ((lj), (|lS|) , and (|l|), we are able to find the 
scalings of surface area and pore space volume with n, 



and 



S(n) 



<&(«) ex n 



,(a-d)/( a -l) 



-/3 cx n _(d+1_Q)/(Q ~ 1) . 



(22) 



(23) 



These relations provide us with further methods for 
determining a and for testing scaling relations between 
the iterative and structural nature of the packing. Equa- 
tion |23|, in particular, affords a more robust measure- 
ment than does obtaining a directly from P(r;n), and 
we employ this fact in our numerical investigations. 



Note that for comparison between the current theo- 
ry and that of ABK we must convert scaling predictions 
that depend on time to those that depend on sphere num- 
ber. Specifically, the ABK theory predicts $(£) oc t~ A 
where A = exp[2 - (2 d+2 - 2)(d + 2)]. Because the 
ABK model can be seen as an RAP model with infi- 
nite growth velocity, the time between events is inversely 
proportional to the pore space available (given a uni- 
form rate for attempted nucleation k). So the result for 
the decay of pore space from ABK may be rewritten as 
<&ABK{n) cx , where 



1 " 1 
-Y - . 

K l^l ®ABK{i) 



(24) 



This implies that t n cx n 1+!3 ' . So if ABK predicts §(t) cx 



r 



this is equivalent to writing <&abk (n) 



-A{l+f}') 



Equating the two forms leads to the conclusion that 

A 



I- A 



(25) 



This will be useful in the following section when the scal- 
ing of pore space decay is examined. Finally, key scaling 
relations are summarized in Table Q. 



relation 


exponent 


P(r) oc r' a 


Q 


r c (n) oc n~ s 


5 = d+ 1 - a 


S(n) cx n 1 


j = (a - d)/(a - 1) 


$(n) cx rT 13 


j3 = (d+ 1 - a) /{a - 1) 



TABLE I: Scaling relations between exponents for various 
parameters valid for general packing-limited growth models: 
P(r) is the probability distribution of radii, r c is the radius 
of the typical largest sphere that may be inserted, S(n) is 
the surface area of n packed spheres, and &(n) is the pore 
space volume. All exponents are given in terms of a and the 
dimension d. 



V. NUMERICAL RESULTS 

The numerical procedure for all variants of PLG 
schemes obey the same basic algorithm. Spheres are seed- 
ed according to a Poisson distribution, with rate k and 
only grow within the pore space. Once injected, spheres 
grow according to model rules (linear velocity, hetero- 
geneous, exponential, etc.) until they collide and are 
jammed. In the case of the RAP model, only one sphere 
is allowed to nucleate and expand at any given point in 
the simulation. Using this event driven procedure we are 
able to calculate limiting states and determine improved 
estimates of the universal value of a for d = 2, 3, and 
4. All packings are simulated on a unit hypercube with 
periodic boundary conditions. Figure [| shows an exam- 



FIG. 3: A random Apollonian packing initially seeded with 
two circles after 10,000 circles have been placed. The density 
of the system is p ~ 0.94. 



pie packing using the RAP model in d = 2. It is clear 
that the number of spheres will increase without bound 
and that the pore space ultimately vanishes. The associ- 
ations with the Apollonian case are visually evident. 

RAP Heterogeneous Exponential Linear 



k n/a 
Gi(t — ti) oo 
N s 50 
JV tot 2.5 x10 s 
p 0.97 



10 

0.2 - 2.0 

48 

2.3 x 10 6 
0.96 



10 

e t-U 

10 

5.4 x 10 5 
0.91 



10~ 5 
0.2 
1 

3.4 x 10 5 
0.97 



TABLE II: Simulation details for four different PLG models. 
k, is the per unit volume sphere nucleation rate, d(t — ti) is 
the growth rate of the the ith sphere which is nucleated at 
time ti, Ns is the number of simulations, N^ Q ^ is the total 
number of spheres placed, and p is the approximate limiting 
density of each simulation. The sphere radii distributions for 
these models are shown in Figure U. 



To demonstrate the universality of a, we consider four 
different PLG models in d = 2 with a variety of initial 
conditions, the details of which are given in Table 0. The 
frequency distribution N(r) for these four PLG models is 
shown in Figure ^. The main plot shows the distributions 
recentered using their respective means. The recentered 
distributions are indistinguishable to the eye clearly indi- 
cating that the specifics of the growth mechanisms are 
irrelevant and that the exponent a is a universal one. 

The scaling theory developed in the previous section 



FIG. 4: The form of N(r) for four models: random Apol- 
lonian packing (circles), heterogeneous linear growth rates 
(diamonds), exponential growth rates (triangles), and linear 
growth rates (squares). Simulation details are contained in 
the text. The inset double-logarithmic plot shows N(r) ver- 
sus r for each model shifted horizontally for the purpose of 
illustration. The slope of the solid straight line indicates the 
d = 2 universal exponent a ~ 2.56 (see Table III). 



relates the structure of a packing (its fractal dimension) 
to the change of basic elements (critical radius, surface 
area, volume, etc.) as a function of iteration. In order to 
test the scaling theory, as well as the specific prediction 
a = d+ (d+l) / (d+2) in higher dimensions, we now exam- 
ine iterative and structural scalings of the RAP model in 
d = 2, 3, and 4. 

We estimate a using a var iety of means and the results 
are summarized in Table III . The agreement between the 
predicted value of a calculated using <I>(n) and that calcu- 
lated directly from the geometry of the packing provides 
further proof that the scaling theory outlined in the pre- 
vious section is valid, regardless of the specific value of a. 
The current theory appears to improve with increasing 
dimension, whereas the approximate theory of ABK only 
holds in d = 2. 

Both the current theory and that of ABK predict the 
scaling of pore space as a function of iteration or time. 
Evaluating the pore space decay $(n) in d — 2 suggests 
that the agreement with ABK is most likely coincidental. 
As evidenced by the data in Fig. |^ we are able to establish 
95% confidence intervals that exclude the predictions of 
ABK in d = 2. This is not altogether surprising as ABK 
is an approximate theory describing the collision of mov- 
ing spheres whereas in the RAP procedure all collisions 
are between a sphere and the presently jammed state. 
In higher dimensions, the theory of ABK fails, and the 
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d = 2 


d = 3 


d = 4 


$(n) 


2.564(1) 


3.733(2) 


4.833(2) 


$(r) 


2.58 


3.73 


4.83 


S(r) 


2.55 


3.74 


4.86 


JV(r) 


2.53 


3.70 


4.79 


a = d+[(d+l)/(d + 2)] 


2.75 


3.8 


4.83333 


ABK theory 


2.553740 


3.94505 


4.99904 



TABLE III: Estimates of a using various numerical methods 
and comparisons to theory. All numerical results are from 
simulations containing 5 x 10 6 spheres. Parentheses indicate 
95% confidence intervals on last digit. 



-0.5 



O 




current theory becomes increasingly appropriate, falling 
within the numerically defined confidence intervals for /3 
in d = 4. The results are summarized in Table [V. 
VI. CONCLUSION 



The iterative nature of the random Apollonian pack- 
ing model and its simple prescription for packing spheres 
suggests that at sufficiently small scales, the structure of 
packing-limited growth models is essentially unaffected 
by initial conditions or dynamics. The scaling of pore 
space, surface area, number, and critical radius r c are all 
interrelated and can be expressed in terms of simple pow- 
er laws using only the dimension d and the exponent a. 
Extensive numeric simulations demonstrate the validity 
of the predicted upper bound a. Although a appears to 
be an overestimate of the actual value of a in d = 2 and 
d = 3, it agrees with simulations in d = 4, and presum- 
ably beyond. A refined approximation of the insertion 
probability P; ns (r;n) seems to be the best approach by 
which to improve the estimate of a. The applicability 
of predictions of the present model and its relevance to 
physical and biological problems may lie in the process of 
balancing the idea of packing-limited growth with other 
dynamical possibilities such as aggregation, competition, 
and death. 



FIG. 5: The decay of pore space volume $(n) in d — 2, 3 
and 4 correspond to the solid, dashed, and dot-dashed lines 
respectively. The fits to the power law decays are summarized 
in Table | 



d 



d- 



d- 



P 

current theory (23|) 
ABK theory dg5) 



0.278(1) 0.0975(2) 
0.1429 0.07143 
0.2872 0.01832 



0.0434(2) 
0.04348 
2.404 x 10" 



TABLE IV: The predicted exponent for the decline of volume 
fraction, <E>(n) tx n _/3 . The numerical estimates of j3 are taken 
from simulations containing 5 x 10 6 spheres. The decay of 
pore space as a function of iteration can be seen in Figure ^[ 
Parenthesis indicate error on estimate of last digit. 
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